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The vortex phase diagram in the external rotation fre- 
quency versus temperature is calculated for dilute Bose- 
Einstein condensed gases. It is determined within the 
Bogoliubov-Popov theory for a finite temperature where the 
condensate and non-condensate fractions are treated in an 
equal footing. The temperature dependences of various ther- 
modynamic instability lines for the vortex nucleation are com- 
puted to construct the phase diagram. Experiments are pro- 
posed to resolve a recent controversy on the vortex creation 
problem associated with the quantized vortex observation in 
87 Rb atom gases. 

PACS numbers: 03.75.Fi, 05.30.Jp, 67.40.Vs 



I. INTRODUCTION 

Since its realization in 1995 Sj , Bose-Einstein con- 
densation phenomena (BEC) in alkali atom vapors (Na, 
Li and Rb) have been attracting much attention. Sev- 
eral outstanding fundamental questions such as the rela- 
tionship between the origin of superfluidity and Bose- 
Einstein condensation can now address directly both 
experimentally and theoretically Recently long- 

sought quantized vortex has been successfully created in 
87 Rb atom gases by two different experimental methods 
]fj|-(io||; One is the so-called phase imprinting ||f7| and 
the other is to use an optical spoon @-|lc[|. The lat- 
ter is conceptually similar to the rotating bucket exper- 
iment in 4 He and 3 He to create quantized vortices [ Q| ; 
Here the harmonically trapped BEC with an elongated 
cigar shape is rotated around the long axis by stirring 
two laser beams. As the stirred rotation angular fre- 
quency f2 increases, a quantized vortex appears at the 
nucleation frequency Q nuc , followed by multiple vortices 
with each carrying the quantized unit circulation. They 
form an Abrikosov-like lattice structure. The resulting 
nucleation frequency f2 nuc relative to the radial confining 
potential lu± is fl nuc ~ 0.68wx which is not much depen- 
dent on their experimental parameters such as the atom 
number, or the confining potential, etc. in their exper- 
iments [p[-^0[. This particular number fl nuc ~ 0.68lu± 
sharply disagrees with the prediction by the ther- 
modynamic critical frequency ^gj b a ; at which the single 
vortex state becomes stable thermodynamically, namely, 
a theory based on Thomas-Fermi approximation gives 
^ 'global I ' u 'J- ~ 0.37. R turns out that the nucleation fre- 
quency happens to coincide with the characteristic fre- 
quency Cl w =o at which the non-vortex state ultimately 



becomes unstable, leaving only the single vortex state 
stable globally in the energy configuration space. How- 
ever, it is argued [jl3 14 that since the actual nucle- 
ation process is more or less related dynamical instability 
rather than the above thermodymanic quantities, that is, 
the process is initiated by resonantly exciting the quadra- 
pole mode in a slightly radially distorted cigar-shaped 
BEC. Thus there is no definite answer to explain Q nuc at 
present. 

The purpose of the present paper is to complete the 
vortex phase diagram in the plane; the external rotation 
frequency f2 versus temperature T in order to help estab- 
lishing vortex properties in general in the present dilute 
Bose gas systems, in particular, the nucleation process by 
identifying the following characteristic critical frequen- 
cies related to various thermodynamic instabilities, which 
will be found to have different T-dependence. 

Previously in our series of papers Jl5| |l9|] , we have iden- 
tified several characteristic critical frequencies; Upon in- 
creasing from zero, the single vortex state first becomes 
locally stable at Qf oca i and then stable globally at Qgi oba i 
relative to the non-vortex state. At further higher fre- 
quency CI w =q the non- vortex state exhibits an instability 
toward the vortex state and finally the single vortex state 
becomes unstable, signalling that multiple vortex state 
can be stabilized. Here we will investigate the tempera- 
ture evolution of these critical frequencies by solving the 
Bogoliubov-de Gennes equation for non-condensates. It 
is coupled with the generalized Gross-Pitaevskii equation 
for the condensates within the so-called Popov approxi- 
mation. The condensate and non-condesates are treated 
in an equal footing. We assume a cylindrical symmetric 
system confined radially by a harmonic potential. The 
rotation occurs around the cylindrical axis. 

Stringari f2(f | gives a phase diagram in the plane; fl ver- 
sus T and determines the T-dependence of Slg; obQ ; in our 
terminology, basing on the Gross-Pitaevskii equation by 
neglecting the non-condensate contribution. In contrast, 
we take into account thermally excited non-condensate 
contribution fully in our calculation as mentioned. 

This is one of our continuous efforts to understand the 
vortex matter in BEC. We have first pointed out the 
presence of the anomalous negative mode with the par- 
ticular angular momentum both for the bucket-like sys- 
tem with a rigid wall boundary [ fl5| and for the system 
with the confining harmonic potential Jl6| . It turned out 
later that this anomalous mode in BEC is crucial to un- 
derstand the vortex matter, in particular, the nucleation 
problem. Since we find that the non-condensate contri- 
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bution due to finite temperature effect helps stabilizing 
the vortex state by pushing up the negative mode, it is of 
interesting to know the vortex phase diagram in £1 versus 
T. 

The arrangement of the paper is as follows: The for- 
mulation of the problem is given within the Hartree-Fock 
Bogoliubov framework in next section. Based on the self- 
consistent calculations the local stability of a vortex is 
considered in III and the global stability is discussed on 
IV. The final section is devoted to constructing the phase 
diagram in fl vs T and discussions. 



II. FORMULATION 

The Hamiltonian in a rotating frame with the angular 
frequency with u> is given by 

H rot = H-u- J "*t(r)-(rxp)-*(r)dr (1) 

r ~ h 2 

H = / dr*t( r ){ y 2 + V(r) - p, }*(r) 

J 2m 

+ |y dr&(r)&(r)$(r)V(r) (2) 

where the creation and annihilation operators of the Bose 
particle are ^(r) and 'J'(r). They interact each other 
with the contact interaction g = 4 ^ - (a: the scattering 
length and m: the mass). The confining potential V(r) 
and p is the chemical potential. These are decomposed 
into ^(r) = (j>(r) + ip(r) in terms of the condensate wave 
function: 0(r) and the non-condensate part: V'(r). By 
this decomposition, the above Hamiltonian is rewritten 
as 

Hrot = Jdr[{ 4>*{v)h(v)4>(v) + l -gW)\ 4 } 

+ ft(r){h(r)+g\ct>(r)\ 2 }<P(r) + h.c. 
+ ^{v){h{v) + 2g\ ( t > {v)\ 2 }^{v) 
+ ^ {v)^ {y)4>{y)4>{y) + h.c. 

+ g^(r)^(r)tp(r)(j)*(r) + h.c. 
+ |^( r )^t (r) ^ (r) ^ (r) 

+ ih{ (f>*(r)uj • (r x V)^(r) 
+^ t (r)w • (r x V)0(r) 
+<f>*(r)ti ■ (r x V)^(r) 
+^( r )dJ.(rx V)^(r)}]. (3) 

The single particle Hamiltonian is 



h(T) 



2m 



V 2 + V(r) - ft. 



(4) 



^(r)^(r)^(r) ~ 2(7/> t (r)^(r)}^(r) 
+ (^(r)^(r))-0 t (r) 
^(r)^(r)4>(r)ij;(r) ~ 4(V> t (r)?A(r))^ t (r)^(r) 
+ (^ t (r)^ t (r))^(r)^(r) 
+ (^(r)^(r))-0 t (r)^ t (r) 
^(r^r^r) ~ 2p(r)^(r) 
^(r)^(r)^(r)^(r) ~ 4p(r)^(r)^(r) 

with the non-condensate den- 

sity defined as: (ijr (r)tp(r)) — p(r). Thus we arrive at 
the following effective one-body Hamiltonian within the 
present mean-field approximation: 

Hrot = Jdr [ { 0*(r)/i(r)^(r) + | 9 |0(r)| 4 } 

+ ^t( r) { h{r) + ff |0( r) |2 + 2gp {v) }cf>(r) + h.c. 
+ ft(r){ h(v) + 2g\^(r)\ 2 + 2gp{v) }^(r) 
+ |^+(r)-0+(r)^(r)^(r) + h.c. 
+ ih{ (f)*(r)uj ■ (r x V)0(r) 

+^(r)u) ■ (r x V)0(r) 

+</>*(r)w • (r x V)-0(r) 

+-0 f (r)cJ • (r x V)-0(r) } ] 



(5) 



This can be diagonalized by the usual Bogoliubov trans- 
formation by noting the commutation relation: 

[$(r>),H rat ] = { h(r)+g\<f ) (r)\ 2 + 2gp(r) 
+ihu) • (r x V) }</>(r) 
+ h(r)i>(r) + 2<#(r) | 2 ^(r) 
+ g<j, 2 {Y)^{v) + 2gp{v)ij{v) 
-ihuj-[Vip(r) xr]. (6) 

The Bogoliubov transformation is given by 

V-(r)=^[ Uq (r) ??q -^(r) ?7 t]. ( 7 ) 



The resulting diagonalized form reads 
H rot = E + £ qVlVq- 



(8) 



Here we employ the following Bogoliubov-Popov approx- 
imation: 



This diagonalization can be accomplished under the fol- 
lowing conditions. First we impose the condition on <j)(r) 
as 

{ h(r) + 5 |0(r)| 2 + 2gp(v) + ihuj-(r x V) }0(r) = (9) 

which is nothing but the so-called Gross-Pitaevskii (GP) 
equation containing the non-condensate contribution. 
We must fulfill the following equations for diagonaliza- 
tion: 
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{ h(r) + 2.g|0(r)| 2 + 2gp(r) } Wq (r) - # 2 (i> q (r) 

- ifiuj ■ {Vu q (r) x r} = e q u q (r) (10) 

{ h(r) + 2.g|0(r)| 2 + 2gp(r) } Wq (r) - 5 0* 2 (r)u q (r) 

+ ihu; • {Vw q (r) x r} = -e q u q (r). (11) 

This coupled equation is known as Bogoliubov-de Gcnncs 
(BdG) equation |2lj used both for BEC and for super- 
conductors. 

We notice that the eigenfunctions; w q (r) and w q (r) sat- 
isfy the otho-normalization requirement or the positive 
norm condition: 



dr[ u* (r)u q (r) - <(i> q (r) ] = 5 



p,q- 



(12) 



The non-condensate density is determined 
self-consistently by 

p(r) = (^(r)^(r)} 

= Mr)| 2 /( £q ) + l« q (r)| 2 {/( £q ) + 1} ]• (13) 
q 

These equations (9), (10) and (11) constitute a complete 
set of the non-linear coupled equations. 



0, ±2ir/L, ±4tt/L, • • ■, where L is the period of the length 
along the z-axis. 

By choosing the external rotation Q = (0,0, a;), we 
obtain the GP equation as 

r Ti 2 d 2 Id w 2 , TT , . 

2m r dr r z 

+ g(f> 2 (r) + 2gp(r) -hu>w}(/)(r) = 0. (18) 

Similarly, the BdG equations are now 



h 2 d 2 + 1 d {q e + w) 2 
2m dr 2 r dr r 2 



+2g((f) 2 +p)- fiu{w + qg) ]u q (r) 
-g4> 2 (r) VcL {r) = e q u q (r) 
. h 2 d 2 Id (qg- w) 2 



(19) 



2m dr 2 r dr r 2 
+2g(<f> 2 +p) + huj(w - qg) ]« q (r) 
-g<j> 2 (r)uq(r) = £ q w q (r). 



(20) 



The actual implimcntation of the numerical computation 
for these coupled equations is described before [ p7| , ^8[ . 



A. Cylindrical system 

In order to describe a single vortex threading through 
the center of the cylindrical symmetric system, it is con- 
venient to use a cylindrical coordinate: r = (r, 6, z). We 
impose the boundary condition that the condensate wave 
function and the eigenfunctions for the non-condensate 
vanish at r = R, namely, 4>(r — R) — 0, u q (r = R) = 0, 
and Vq(r = R) = 0, where R is taken far enough from 
the center. 

The confining potential is given by 



V(r) = ~m{2m> r ) 2 r' 2 



(14) 



and the periodic boundary condition is employed along 
the z direction. The condensate wave function is written 
in the form: 



(r,9,z)=<t>(r)e iwe . 



(15) 



The winding number w is finite when a vortex is present 
and vanishes when a vortex is absent. In the following 
we consider the vortex with w = 1 only. 

The eigenfunctions are also written by the symmetry 
reasons as 



Uq(r) = u q (r)e i9 « V ( " +tu)e 
w (r) = v^{ry q * z e l{qe ~ w)B . 



(16) 
(17) 



The eigenvalues q = (q r ,qg,q z ) may take the following 
values: q r = 1, 2, 3, , • • -, qg — 0, ±1, ±2, , • • •, and q z = 



B. Free energy 

The internal energy E of the system is calculated by 
taking the Hamiltonian within the mean-field approxi- 
mation, namely, 

E = (H rot ) + fj,N 

= J dr[{0*(r)Mr)<Xr) + ^(r)| 4 } 

+ (^(r){ h(r) + <?|0(r)| 2 + 2gp(r) }<£(r)) + h.c. 
+(ft(r){ h(r) + 2 5 |</>(r)| 2 + 2gp(r) }^(r)> 
+|(^ + (r)^t(r)0(r)0(r))+/i.c. 

+ih{ {<f>*(r)i3 ■ (r x V)</>(r)) + {i)\r)Q ■ (r x V)0(r)} 
+(4>*(r)C3 • (r x V)V-(r)) + (ft(r)<3 • (r x V)^(r)) } ] 
+pN (21) 

or, it can be written as follows: 

E = J dr[{^)ft(r)#r) + ^(r)| 4 } 

+#(r){ h(r) + 2 5 |0(r)| 2 + 2.gp(r) }^(r)} 
+|(^t(r)^t( P )^( r )^ r )) +h . c . 



+i%{4>*{r)u ■ (r x V)0(r) 

+ (^(r)uj ■ (r x V)^(r))} ] + pN 

E + pN + 0!) + pN nc 



(22) 



where we define 
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Eo = J dr [0*(r)/.(r)0(r) + ig|0(r)| 4 

+i%(j)*{r)uj- (r X V)^(r) ] (23) 
H x = /dr [ ^t( r ){ ft ( r ) + 2 3 |0(r)| 2 + 2gp(r) }^(r) 



+ |^( r )^t (r)0(r)0(r) + |^( r )^(r)^( r )0*( r ) 



-M7iV^(r)£J • (r x V)-0(r) 



(24) 



In terms of the cylindrical coordinate r = (r,6,z), we 
obtain 



The area density along the z axis can be estimated as 
n z = 2 x 10 9 /m ~ 3 x 10 9 /m at the center of r = as 
seen from Fig. 1(b). In the following the dimensionless 
"gas parameter" an z = 11 is chosen in the present two- 
dimensional disk system and the rotation frequency is 
normalized by the radial confining potential lu± = 2irv r , 
denoted by SI = ^ from now on (u>±/2ir = llbHz). The 
system size is set to R — 7fim and L = 5/im. The tran- 
sition temperature of the present stsem T c q = 207.8nK 
at O = 0. The energy is scaled by Hlo± = hv r . 



E a = J dr[ 0*(r) { -— V 2 + V(r) - fi } <f>{r) 

+ \gW)t + -Tiujw4> 2 {r) ] 
f , r , * / \ / Ti 2 . r d 2 Id w 2 . , , . 

+V(r)cj) 2 (r) - n<j> 2 (r) + -g^{r) - hujw(j) 2 (r) ]. (25) 

As shown in Appendix, the total internal energy is now 
evaluated by the formula: 

E = (H rot ) + flN 

-™/*Wr)(-|;H £ + WW 

+V{r)<p 2 (r) + -g^ 4 {r) - hujw<p 2 {r) } 

+ E £ q/( £ q) + E e i / dr K\ 2 + v I <M r ) (26) 

q q 

where /(e) = l/{exp(e q /A;BT) — 1} is the Bose distribu- 
tion function. 

The entropy S of the system is obtained by the follow- 
ing expression ||^] : 

S = k B J2i /feq/(e q ) - H 1 - cxp(- £q /fc B T) ] }. (27) 



From eqs. (26) and (27), we can estimate the free energy 
F = E — TS of the vortex and non- vortex states. 



C. Calculated system 

In order to perform the self-consistent calculation to 
determine various characteristic frequencies as a function 
T, we must specify the parameters: As a typical exam- 
ple, we take up the vortex nucleation experiment on 87 Rb 
atom vapor by Chevy et al M. The scattering length 
a = 5.5 x 10~ 9 m and atom mass m = 1.44 x 10~ 25 kg 
for Rb atom. The radial and z axis trapping frequen- 
cies are 175Hz and 10.3Hz respectively. We evaluate the 
density profile of the condensate in Fig.l by solving the 
GP equation at T—0 for the total number N — 2.5 x 10 5 . 



III. LOCAL STABILITY 



A. Local stability of the vortex state 

The local stability of the single vortex state is de- 
fined by the frequency region bound by the two frequen- 
cies: Qf ocal < SI < QY oca [ where the excitation spec- 
trum is well-defined and no negative frequency exists. 
We show in Fig. 2 the density profiles for the condensate 
(a) and non-condensate (b) and also the associated exci- 
tation spectrum (c) for various SI values under a fixed T 
(T/T c o = 0.4). As S7 increases, we can see from Fig. 2(a) 
that the condensate is depressed from the outside region 
while the non-condensate increases its density at the vor- 
tex center, filling the vortex core by the non-condensate, 
eventually a large amount of the non-condensate becomes 
piling up the outside region for larger Q values. The 
filling-in of the non-condensate effectively helps pushing 
up the so-called anomalous mode, that is, the negative 
eigenvalue at qg = — 1 toward a positive one as seem from 
Fig. 2(c), resulting in the local stability of the vortex state 
and giving rise to S7^ ca/ [ p3| . In Fig. 2(a) as the inset 
we display the trace of the negative mode and also the 
lowest mode at qe > which becomes negative at Sl^ ca/ 
at higher Si's. It is seen that for Slf oca i < SI < Sl^ ca/ the 
excitation spectrum is well defined, implyig that a vortex 
is locally stable in the energy configuration space. 

The piling-up of the non-condensate at the outside re- 
gion arises from the repulsive interaction between the 
condensate and non-condensate parts. It is rather re- 
markable to notice from Fig. 2 (a) in this connection that 
the vortex core size measured by the peak position of the 
condensate density is hardly changed under varying SI. 

As T is varied under a fixed SI, both density profiles 
of the condensate and non-condensate change into each 
other shown in Figs. 3 (a) and (b) respectively. As seen 
from Fig. 3(a) the core radius defined by the maximum 
position of the density profile of the condensate slightly 
increases as T increases. The non-condensate density 
quickly fills in the central core region and also the outside 
region as T increases as seen from Fig. 3(b). Again the 
the negative eigenmode at qg = —1 becomes positive at 
^focai' shown in Fig. 3(b) as the inset. It will be seen 
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shortly that Qf ocal is a monotonically decreasing function 
of T. 

When we perform the self-consistent calculation, we 
exclude modes with negative energy in eq.(13). The neg- 
ative mode(s) represent the instability of system consist 
of condensate and non-condensate. The steep step in the 
insets of Fig. 2(a) and Fig. 3(b) appears when the lowest 
mode becomes ignored. 



B. Local stability of the non- vortex state 

The corresponding local stability can be examined for 
the non-vortex state in a similar way: We show the cor- 
responding example for that case in Fig. 4 where (a) and 
(b) exhibit the density profiles for the condensate and 
non-condensate respectively. Again as Q increases, the 
former (latter) is depleted (piled up) from the outside re- 
gion. There is no filling-in effect of the non-condensate 
in this vortex-free state. 

The excitation spectrum as depicted in Fig. 4(c) shows 
the appearance of the negative mode at 0^=0 for certain 
mode with qg > 0. This signals an intrinsic instability 
of this non-vortex state toward some other state, most 
likely toward the vortex state in the present situation. 
The trace of the lowest mode with qg > is shown in the 
inset of Fig. 4(a) which is seen to Q w -o = 0.566 in this 
case. 

The T-dependent poperties of this vortex-free state are 
depicted in Fig. 5: As T increases, the density conversion 
from the condensate to the non-condensate occurs. This 
is shown in Figs. 5(a) and (b) where each density dis- 
tribution is displayed. The inset of Fig. 5(a) exhibits the 
trace of the lowest mode with qg > 0, indicating that the 
negative mode in qg > become positive at T/T c q = 0.4 
in this example. 



IV. GLOBAL STABILITY 



with no the peak density, R the trapping radius and £o = 
l/87rano. The angular momentum L per unit length is 
given by L = UN, thus the critical frequency Qg lobal is 
determined by 



^globalL — 0- 



(30) 



Under the Thomas-Fermi approximation this leads to 



^ global ~ Z mR 2 



In 



0.888.R 



(31) 



which indicates the global stability of a vortex relative to 
the non- vortex state, where R is the radius determined by 
the Tomas-Fermi approximation. At T=0 we can easily 
estimate fogi b a i — 0.372 for the Thomas-Fermi approxi- 
mation, which is in good agreement with our own calcu- 
lation Qgi oba i = 0.382. Stringari (2(J further extends it 
to finite temperatures by neglecting the non-condensate 
contribution. We have performed the fully self-consistent 
calculation, taking into account the non-condensate con- 
tribution. 

We compare various approximate estimates Qgi oba i as 
a function of T. First the Stringari's ^ig [oba i is well re- 
produced by our numerical calculation using the expres- 
sion (31) where no(T) and R{T) are calculated by our 
self-consistent solution. Starting with eq.(14) and equat- 
ing the corresponding non- vortex energy, we can evaluate 



global 



as shown in Fig. 6 where curve (2) indicates the 
result when neglecting the non-condensate contribution 
while (3) is fully consistent results. It is seen that as T in- 
creases two curves deviate progressively showing the im- 
portance of the non-condensate contribution. However, 
curve (3) neglects the effects of the entropy S defined by 
eq.(27). In fact, the global stability of the vortex state 
relative to the non-vortex state at finit temperature is 
defined by the comparison of the two free energies in the 
system with and without a vortex. The T-dependence of 
the critical value is plotted as curve (4) in Fig. 6. 



The relative energies of the non- vortex and the vortex 
states should change as functions of and T. Lundh et al 
|jl2"f estimate the energy increment e per unit length due 
to the vortex creation for the uniform system at T = 
as 



Ti , 1.4646 
e = Trn — In — - — 

m 4 



(28) 



where n is the condensate density, b is the radius of a 
bucket and the coherent length £ = l/87ran. In the pres- 
ence of the confining potential the above expression is 
modified to 



Tf , 0.888i? 

ttuq — In — 

m £ 



(29) 



V. PHASE DIAGRAM AND REMARKS 



We show in Fig. 7 the four characteristic frequencies; 
^Laii n global, n ^=o and nf ocal together with the critical 
temperature T c (fi). This is estimated as 



T c (fi) = (1 - tt 2 )3 



(32) 



by Stringari p0| for a non-interacting system, which 
is also approximately valid for the present dilute Bose 
gas. It is seen that ^i oca i is a decreasing function of 
T because the non-condensate fraction excited thermally 
serves as stabilizing a vortex, a fact pointed out already 
by Isoshima and Machida fig ], and also recently by Vir- 
tanen et al p3[ . 0^ obai (T) is almost T-independent up 
to 0.5T c o and then quickly increases toward the T c (f2) 
line. According to eq.(31) they meet tangentially. 
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As we have discussed on III, the instability of the non- 
vortex and the single vortex state is defined by the same 
reason. However, as seen from Fig. 7, the critical value 
of the instability the critical Q w =0i which correspond to 
the instability of the non- vortex state, disagrees with the 
critical SlfocaV We interpret the disagreement by the fol- 
lowing effect. From the eigenvalue equations (19) and 
(20), we define the effective potential 

V ef f(r) = V(r)+2g[\ ( f ) (r)\ 2 +p(r)}. (33) 

In comparison with the non- vortex state, the density dis- 
tribution of the condensate |0(r)| 2 in the vortex state is 
slightly pushed out because of the angular momentum. 
As seen from Fig. 8(a), this situation is fully understood 
in terms of V e ff(r). As discussed in III, the eigenval- 
ues with qe > move down as Q increases. This means 
that the non-condensate quickly grows at the edge of the 
condensate, and the non-vortex state (the single vortex 
state) becomes unstable at Q^o (^Yocai)- In short, the 
higher requires in order for the non-condensate to in- 
crease at the edge of the condensate in the vortex state, 
and we find that SlYocai ls always higher than Q w= q. 

The remaining SI w= q and fi^ j are a rather strong 
increasing function of T. As mentioned before, £l w= o 
turns out to be very near the observed nucleation fre- 
quency by Chevy et al [||, which is independent of the 
details of the experimental parameters used. They find 
Qnuc ~ 0.68wj_. This is compared with our Cl w= o = 0.55, 
which is nearer than the thermodynamic critical fre- 
quency Qg [oba i ~ 0.39 at T = 0. This may be pure acci- 
dent, but we can check it by quantitatively examining it 
through its T dependence. 

It is important to notice that ^Y oca i(T) is not much 
different from £l w =o(T) in its value as seen from Fig. 7. 
Thus once a vortex nucleates at Q w= q(T), slight increase 
of SI brings multiple vortices at ^ ca ;(T), that is, the fre- 
quency window of the single vortex state is rather narrow. 
This is indeed an experimental fact || (Cl nuc /u)x — 0.668, 

S^multi— vortices j ^ _L 0.697). 

If the system is kept rotating during cooling from above 
T c , then the nucleation likely begins when Q reaches 
Slf ocal . If the system, cooled already, on the other hand, 
begins rotating, then at fl w =o a vortex may appear. It is 
known for some time that in superfluid 4 He the former 
is lower than the latter as expected while in BEC both 
nucleation protocols give rise to almost same nucleation 
frequency. This is quite understandable because two fre- 
quencies; Slf oca i and Sl w= o become almost equal at T = 
when the interaction is weak while these two frequencies 
become quite distinctive when the interaction is strong 
(see Fig. 3 in Isoshima and Machida p8|). It is also true 
for finite temperatures that the two characteristic fre- 
quencies move oppositely upon increasing T as seen from 
Fig. 7. Here there is a good chance to directly confirm 
the characteristic nucleation frequencies for two different 



cooling protocols, a situation quite different from the su- 
perfluid He rotating bucket case where the rough surface 
of a bucket prevents us from knowing the intrinsic nucle- 
ation frequency. The present gaseous system is free from 
this problem. 

Quite recently, yet another vortex experiment is re- 
ported by Abo-Shaeer et al p3] where many vortices 
over 100 are observed. The very low nucleation frequency 
Q nuc ~ 0.25wx in their experiment is qualitatively un- 
derstandable because the gas parameter an z is large (an 
order of 500), implying that the characteristic frequen- 
cies Sl w= o, Qgiobai and Slf ocal are all small in that region 
(see again Fig. 3 in Isoshima and Machida fTsf] ). 
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APPENDIX: 

We first express each term of eq. (4) in terms of rf 
and r\ as 

^(r){ h(v) + 2g\<p(r)\ 2 + 2gp(v) }^(r) 

= Et ifa i <( h + 2 #H 2 + 2 9P> 3 } 

+mVj {v l {h + 2g\4>\ 2 + 2gp)v* } 
-4r,\ {u*(h + 2g\cj ) \ 2 + 2gp)v* } 
-mm {<{h + 2g\<t>\ 2 + 2gp)u J ■}] 
ft(r)ft(r) = Y^i vlv}u*u* + THVjViVj - 2i 1 \r lj u*v ] ] 

4>(r)4>(r) = El ViVjUiUj + r]Wj v i v j ~ ^ViVjUiV* } 
i,j 

ih { $t(r)cD • (r x V)^(r) } 

= ih^2(u*r]l - Vir)i)u ■ (r x V)(u,%- - v*g]) 

= -ih^2(u*r]} - v^uj ■ [ V(Ujr)j - v*^) x r ] 

= -ihj^uj- [ r)iVjut(Vuj x r) - g\rj)u*{yv* x r) 

-rjirjjViCVuj x r) + TftT/ju^Vu- x r). 
We substitute these terms into eq.(4) and obtain 
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Hi = I dv [ vh { <{h + 2g\<\>\ 2 + 2gp) Uj 
id J 

-Hm] {v*{h + 2g\4>\ 2 + 2gp) Vj 



-~UjVi(f> 2 - ^UiVj^* 2 - ihv*uj ■ (Vvj x r) } 



+(r)iT]jterm) + (rjjmterm) 



(Al) 



where the coefficient of the 77^77^- 

9 » 

-ihu*u ■ (Vuj x r) (A2) 



u*{h + 2g\(j)\ 2 + 2gp)u j --u*v j 4> 2 - -ujv*(f>* 2 



is rewritten by using 

{ h + 2g\0\ 2 + 2gp } Uj - g\<P\ 2 Vj 

—ihuj ■ (Vuj x r) = EjUj 



(A3) 



as 



i2„,*„. 9 * ±2 9 * 1*2 



(8) = EjUjUj + £/^> ~ ^Vj(p' - -UjV- 

= £ju>i- (A4) 
Here we have utilized a relation: 

E, * j.2 5 *i*2\ 

= -f «i(0»i(r)^(r) - § «i(r) Bj *(r)^(r) } 



^2(-9<P 2 u*Vj). 



(A5) 



Similarly, the coefficient of the term r\ry in eq.(a7) is 
found to be 



Thus the expectation value of H\ is finally evaluated as 
follows: 



(Hi) = Zy dr t £ j u i u i(vhj) ~ £jvtvj(r)iVj) 

= H £t j dr t M 2 -N 2 -Z £i / dr l^l 2 

= E £ *v7k->-£^ / drM 2 (A6) 



with 



pN nc = p J dr(^t( r )^( r )) 
= M / *p(r). 



(A7) 
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Fig. 1. The solution of Gross-Pitaevskii equation at 
v r = 175-ffz, v z = 10.3Hz. (a) Stereographic view of 
the density distribution of the condensate |</>(r)| 2 . (b) 
The area density n z per unit length along z-axis at the 
center r = 0. The z-direction of the cylindrical container 
is assumed to be uniform with the average area density, 
n z = 2.0 x 10 3 /fim in the following calculation. 



1 5 
„ 1 

i 05 






— 


- — 





1.0 2.0 3.0 4.0 5.0 6.0 
[■ [ftm ] 



(c) 

Fig. 2. At T/T cQ = 0.24, (a) the ft dependence of the 
condensate density |0(r)| 2 . The solid, dashed, dotted, 
dot dashed lines correspond to Q = 0, 0.25, 0.6, 0.7 re- 
spectively, (b) The dependence of the non-condensate 
density p(r). (c) The lowest edge of the eigenvalues along 
qg for = 0, 0.25, 0.6, 0.7. The eigenvalues at qg = — 1 
move up while the eigenvalues in qg > move down as O 
increases. The inset in (a) shows the trace of the eigen- 
values at qg = — 1 and qg = 7 as a function of f2. 
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(b) 

Fig. 3. At O = 0.2, the T-dependence of (a) the 
condensate density \4>{r)\ 2 and (b) the non-condensate 
density p(r). The inset of (b) shows the trace of the 
eigenvalue at qg = — 1 as a function of T. 
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Fig. 4. At T/T c0 = 0.10, the ^-dependence of (a) the 
condensate density |0(r)| 2 and (b) the non-condensate 
density p(r) for the non- vortex state, (c) The lowest edge 
of the eigenvalues along qg for the selected ft = 0, 0.2, 
0.4, 0.6 respectively The positive eigenvalues in qg > 
become negative as Vl increases. The inset of (a) shows 
the trace of the eigenvalue at qg = 6 as a function of Cl. 
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(b) 

Fig. 5. At = 0.6, (a) the T-dependence of (a) the 
condensate density ^(j)! 2 and (b) the non-condensate 
density p{r) for T/T c0 =0.05, 0.14, 0.24 and 0.48 in the 
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non-vortex state. The inset in (a) shows the trace of the 
eigenvalue at qg = 6 as a function of T. 
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Fig. 6. Curve (1) corresponds to the critical Qgi oba i(T) 
based on Thomas-Fermi approximation. Curve (2) and 
(3) express the critical Qgi oba i(T) in which the energy 
of the vortex state is equal to the non- vortex state under 
same condition. Curve (2) corresponds to Qgi oba i(T) esti- 
mated by neglecting the energy Eq of the non-condensate 
while curve (3) corresponds to the critical angular veloc- 
ity by comparison the total energy E of the vortex and 
non- vortex state. Curve (4) corresponds to Qgi oba i{T) 
by comparing the free energies F of the vortex and non- 
vortex state. 
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Fig. 7. Phase diagram: fl vs T/T c q. (a) the four angu- 
lar velocities, Q.f ocal , Sl% lobal , fl w=0 , and fi£ 0O j together 
with the critical temperature T C (Q) are plotted as a func- 
tion of T. These four lines divide the whole area into five 
regions, (I)-(V). Each of the regions is explained in (b). 
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(b) 

Fig. 8. At T/T c0 = 0.24, fl = 0.6, (a) the effective 
potential V e //(r) denned by eq.(33) for the vortex state 
(w=l) and the non-vortex state (w=0) and (b) the con- 
densate density |</>( r )| 2 in the vortex state and non- vortex 
state. 
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